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ABSTRACT 

3 3 . 1^3 

This report is concerned with the problem of elastic-plastic bending of 
circular plates with axi-symmetric loading and support conditions. Incremental 
theory of plasticity with von Mises yield condition, and associated flow rule 
together with the Kirchhoffean small deformation theory of plates have been 
adopted. Finite element type of approach using direct stiffness method of 
matrix analysis of structures was employed. An annular plate is taken as the 
primary element which is subdivided along its depth into a number of layers. 

The material properties are assigned to each layer, and stiffness matrix of 
the element has been established. Loadings are applied in increments. 

Examples of simply supported and clamped circular plates made of elastic- 
perfectly plastic as well as elastic-plastic with isotropic strain hardening 
materials are given. Other support conditions and material properties can 
be treated. The effect of change of material properties within each increment 
of loading, number of elements, and number of layers have been studied. 
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INTRODUCTION 

The application of the incremental (flow) theory of plasticity in the 
solution of structural problems , particularly those related to plates and 
shells, has caused many mathematical difficulties. The closed form solutions 
have been obtained for a very few simple problems such as thick wall cylinder, 
wedge, and torsional shaft made of nonhardening materials with one or at 
most two parameter loading. Further restrictions, such as incompressibility, 
have been occasionally imposed to simplify the solution. No general approach 
to the solution of this class of problems exists so far. Each solution re- 
quires a special mathematical treatment and researcher’s ingenuity. 

A significant amount of literature is devoted to the development of the 
constitutive laws of plasticity both in micro and macro scales. The foundation 
of the theory is not as yet firmly established and further research will be 
needed in the future. 

However, due to the great demand in the design of lightweight structures 
to achieve economical designs and the most realisitic appraisal of the 
behavior, the best available theories have been employed in the analysis. 
Several approaches have been pursued. 

During the last decade considerable attention has been given to the so- 
called limit analysis approach. The bounds of the collapse load, in the sense 
of small deflection theory, have been obtained for a number of problems. 

* See for example: Prager and Hodge, "Theory of Perfectly Plastic Solids," 

John Wiley, N. Y . , 1951. 


3 


Lower and upper bound theorems were employed in attaining the collapse load. 

In order to apply these theorems, the materials should exhibit considerable 
deformation under constant load such as in the case of mild steel. Although 
limit analysis turned out to be fruitful in the analysis of rigidly connected 
frames where strength is the controlling factor in the design, it generally 
does not lead to any physically significant results in plates and shells. This 
is essentially so because of a favorable change in geometrical configuration 
during loading which brings the membrane forces into action. 

A number of solutions are also available for large deformation, for which 
the elastic deformation is quite negligible compared to the unrecoverable 
component of deformation. These are particularly applicable in the metal 
processing. 

The attempt to solve the problems, utilizing the total or deformation 
theory of plasticity has not yet been completely abandoned. In spite of its 
mathematical inconsistency and physical unsoundness, in general, it can be 
used satisfactorily if the stress path in the process of loading can be 
pred i c ted bef orehand . 

For those problems where the plastic deformation sets in with a small 
change in geometry a more realistic solution which takes elastic as well as 
plastic deformation into account should be sought. This will enable a better 
prediction of deflection and stress distribution than the one using elastic 
theory . 
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At the present time, the use of digital computer in the solution of plas- 
ticity problems appears to be very fruitful. It is quite feasible to antici- 
pate that in near future many important problems will be solved with the aid 
of a computer. The improvement of numerical techniques reduces the numerical 
errors and thus increases the reliability of the results. 

This report is concerned with the analysis of circular plates under arbit- 
rary axisymmetric loading. The support conditions are axially symmetric, but, 
otherwise, they can be specified arbitrarily. The examples are given for 
clamped and simply supported solid plates. By the proposed approach, the 
problem of overhanging or annular plates also can be handled without much 
difficulty. The plate may be connected with other structures such as circular 
cylinders or rotational shells, and the developed method can be used to 
determine the redistribution of stresses at junctures due to plastic deformation. 
The material properties used in the proposed analysis can be completely arbit- 
rary. This can be specified from experimental data or any valid analytical 
expressions. To illustrate the method, elas tic-perf ectly plastic and elastic- 
plastic isotropic strain hardening materials were employed. Any other laws or 
experimental data can be specified as desired. 

The procedures discussed here make use of incremental law which enables 
one to trace the loading history. The internal forces and deflections can 
be determined at any stage of loading. The solution is restricted to small 
deformation so that deflection remains small in comparison to the thickness. 

In other words, the influence of membrane stresses which may be developed as 
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the result of deformation is not accounted for. Shear distortion is also 
neglected in conformity with the conventional Kirchhoff 's hypothesis which 
asserts that a straight fiber perpendicular to the middle plane remains 
straight, unstretched, and perpendicular to the deflected state of the middle 
surface after deformation. 

The plate is divided into a number of annular elements. Each element is 
further sub-divided into a number of layers along its depth. The material 
properties were assigned to each layer at every stage of loading. Direct 
stiffness method of matrix analysis of structures was used. 

The effect of the change of number of elements and layers as well as the 
influence of variation of material properties within each increment of loading 
were studied. 

All computations were carried out by IBM 7090-7094 digital computer 
available at the Computer Center of the University of California, Berkeley, 
using FORTRAN IV language. 
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I. REVIEW OF PREVIOUS WORK 


The problem of behavior of axi-symme trica lly loaded and supported circular 
plates loaded beyond elastic limit has attracted the attention of a number of 
investigators, because of the apparent simplicity of the problem and the vast 
number of applications. A brief review of some of the work; utilizing 
deformation theory of plasticity, limit analysis, and incremental theory of 
plasticity together with the results of some of the experimental investigations 
which have been published so far follows. 

1.1 Plate Bending Using Deformation Theory of Plasticity 

The Heneky type deformation theory of plasticity together with Huber- 
Mises yield condition was used by Sokolovsky 1944 [1,2] in the solution of 
bending of circular plates. Kirchhoff’s hypothesis and small deflection 
theory were employed. The problem of simply supported plate under uniform 
and axisymmetrical partial uniform load which contains as its extreme case 
concentrated load at the center was solved. The material is considered to be 
elastic-perf ectly plastic and the scheme to solve for material with strain 
hardening was also indicated. 

Bending of circular and annular plates with variable thickness have 
been discussed by Grigoriev [3], Dvorak 1959 [4] discussed an annular plate 
which is subjected to a ring load at outer boundary and has a simply supported 


inner boundary. 
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The problem of circular plates clamped around the boundary and subjected 
to uniform and partial circular uniform load was solved by Ohashi and Murakami 
[5] 1964. The material is considered to be non-hardening and obeys von Mises ' 
yield condition. 

1.2 Limit Analysis 

The extension of the theorems of limit analysis* encouraged many authors 
to apply them to the problems of plates and shells. 

The early attempt to apply limit analysis theorems to obtain the collapse 
load of circular plates was made by Pell and Prager in 1951 [6]. The materials 
were considered to obey von Mises' yield condition. A simply supported plate 
subjected to uniformly distributed load was discussed, the bounds of the collapse 
load were established, and an approximate value for the load-carrying capacity 
was suggested. Hopkins and Prager in 1953 [7] discussed the problem for 
materials obeying Tresca's yield condition. Simply supported and clamped 
plates subjected to certain simple types of axially symmetric loading were 
treated. The exact collapse loads were obtained. Drucker and Hopkins in 1954 
[8] extended the work to the case of large and small overhang. Hopkins and 
Wang in 1955 [9] compared the collapse loads obtained by utilizing Mises, 

Tresca, and parabolic yield conditions for both simply and built-in supports. 

It is interesting to note that the ultimate load derived by Sokolovsky in 

1944 [l] using deformation theory of plasticity agrees closely with the collapse 

♦Drucker, D.C. , Prager, W. , and Greenberg, H.J., "Extended Limit Design 
Theorems tor Continuous Media." Quart. Applied Math. Vol. 3, No. 4, Jan. 

1952, pp. 381-389. 
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load obtained using limit analysis theorems. Finally, a set of charts for the 
design of circular plates under axially symmetric loading was compiled by Hu 
in 1960 [11]. 

The generalization of the collapse load of circular plate under con- 
centrated load for plates of arbitrary shape was made by Schumann in 1958 [12] 
and also for variable fixity by Zaid in 1958 [13]. It was found that for 
materials obeying Tresca’s yield condition (regardless of shape and end fixity) 
the limit load is equal to 2n times the unit yield moment. 

The load carrying capacities of annular plates with either fixed or 
simply supported edge conditions were discussed by Chernina in 1958 [14], 

Treses' s yield condition was utilized. The case of simply supported edge 
conditions was later corrected by Hodge in 1959 [ 15 ] . 

Attention was also directed towards the minimum weight design of circular 
plates, utilizing limit analysis method. References [16] through [20] can be 
cited here. In 1964 an analog model was devised by Marcal and Prager [21] 
and applied to a circular plate problem. 

Limit analysis was also employed to obtain the collapse load of cir- 
cular plates made of initially anisotropic materials. In 1956 Sawczuk [22] 
discussed a solid plate with either simple or fixed boundary subjected to 
uniform load. The materials of plates were considered to be cy lindrically 
orthotropic, obeying the modified Tresca's yield condition which has 
different platic moduli in two perpendicular directions. Non-homogenei ty has 
also been d iscussed in this paper. Hu 1953 [23] studied the similar problem; 
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and a series of design charts were constructed for ultimate design of circular 
plates under axisymmetric loading by Markowitz and Hu in 1964 [24]. Mura et 
al (1964) [25] extended the work for materials obeying the yield criterion 
suggested by Hill, which reduces to the Mises yield criterion for isotropic 
materials. 

The problem of interaction of in plane tension and bending in annular 
plates was discussed in 1960 by Hodge et al [26], and bounds on the interaction 
curve were established for various inner edge support conditions. 

In 1963 Sawczuk et al [27] studied the effect of transverse shear on 
the plastic bending of simply supported circular plates. Kirchhoff’s hypothesis 
was modified to allow for transverse shear deformation. Both Tresca’s and 
Mises' yield conditions and their associated flow rules were employed. Collapse 
loads were calculated for plates under uniform load over central area. 

Finally, linear programming was also used to obtain the collapse load 
of plate problems [28]. Examples of simply supported and clamped circular 
plates under uniform load are given. 

Detailed discussions of some of the above solutions can be found in books 
by Hodge, 1961 [29], and 1963 [30] , and Sawczuk and Jaeger, 1963, [31 ]. 

1.3 Rigid-Work Hardening Materials 

The use of piece-wise linear yield conditions and associated flow rules 
were suggested by Prager 1955 [34] and Hodge* for the solution of work- 

*Hodge . Jr.. P.G.. The Theory of Piecewise Linear Isotropic Plasticity," 
Deformation and Flow of Solids, Colloq. Madrid, Sept. 1955, Edited by 
R. Grammel, Springer. 
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hardening problems. It allows total stress-strain laws to be used in the small, 
at the same time retaining the characteristic features of incremental laws in 
the large. As indicated by Hodge [35], the plastic flow rules can be explicitly 
integrated under restrictive conditions, defined as a "regular progression. " 

That is, the stress point may not move from one side to another, from one 
corner to a side, or back into the elastic zone. This imposes a serious re- 
striction which does not allow it to be used for problems such as clamped 
plates . 

As an example Prager [34] presented the analysis of a simply supported 
circular plate subjected to a uniformly distributed transverse load. Tresca 
yield condition together with the linear kinematic hardening was used. The 
same problem was discussed by Boyce 1956 [32] using a piecewise yield condition 
which approximates that of Mises f . Later he extended the solution to a 
partially clamped plate [36]. In 1955 Eason [33] discussed the problem of 
plates under a concentrated load. 

The concept of linear, isotropic strain hardening was employed by Hodge 
(1957) [35] in the solution of an infinite annular plate. Tresca yield 
condition was employed and the annulus was assumed to be simply supported at 
its inner edge and free at infinity. It was subjected to a slowly increasing 
moment applied to its inner edge. The loading process was divided into four 
stages, of which the fourth stage deviates from the regular progression. 
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Perrone and Hodge (1959) [37] derived two sets of strain hardening flow 
laws based on kinematic hardening for circular plates made of an initially 
Tresca material. These laws, called complete and direct hardening, differ in 
the point where the plane stress assumption is introduced. The solutions of 
finite as well as infinite annular plates and simply supported circular plates, 
using the derived laws, are reported. The results were compared for complete, 
direct, and isotropic hardening. 

Hwang (1959) [38] treated the problem of simply supported circular 
plates under uniform load. Mises yield condition and associated flow rule 
together with the isotropic strain hardening were employed. Numerical inte- 
gration was used to solve the set of non-linear differential equations. 

Finally, in 1963, Chzhu-Khua [39] using Tresca yield condition and 
linear strain hardening , presented the solution of plate under partially 
uniform load. 

\ ^ Elastic-Plastic Analysis Utilizing Incremental Theory of Plasticity 

An early attempt to estimate the deflections of circular plates made 
of elastic-plastic materials using incremental theory af plasticity was 
made by Hay thorn thwaite in 1954 [40]. The yield condition of Tresca and 
associated flow rule were employed. The key assumption was made that at any 
point on the plate the entire thickness was either fully elastic or fully 
plastic. An annular plate simply supported at the outer edge and clamped 
to a centrally loaded rigid disc at the inner edge was analyzed. The 
results obtained were compared with experiment. 
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Gaydon et al (1956) [41] discussed the problem of circular plate under 
a uniform moment around its edge. Mises yield condition and its associated 
flow rule were assumed to be valid. 

Olszak et al (1957) [42], 1958 [44], [45] treated the elastic-plastic 
bending of non-homogeneous orthotropic circular plates. The material was 
assumed to exhibit no strain hardening and to obey the generalized Mises* 
yield condition. The general moment-curvature relations were derived by 
ordinary plate theory. A restriction that the principal curvatures progress 
proportionally was later introduced to make the expressions more tractable. 
Limitation is similar to that of "proportional loading." The authors report 
that if we require the continuity of stress at the elastic-plastic boundary, 
it is not possible, in general, to satisfy the normality rule which is a 
consequence of Drucker*s postulate of stability. On the other hand, utilizing 
the associated flow rule, it is not possible, in general, to attain the 
continuity of stresses. 

The analysis of a clamped circular plate made of incompressible elastic- 
perfectly plastic materials obeying Tresca*s yield condition and the associated 
flow rule was discussed in 1957 by Tekinalp [43]. An assumption that any 
plate element is either entirely elastic or entirely plastic was made. This 
is strictly valid for a sandwich plate. 

Eason (1961) [46] discussed the problem of simply supported plate under 
partial uniform load. The von Mises yield condition and the corresponding 


flow rule were assumed for the elas tic-perf ectly plastic material of the plate. 
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The stress field was obtained, but the velocity field was derived for the 
limiting case of concentrated load. A comparison has been made with solutions 
obtained utilizing the Tresca yield condition. It is reported that the stress 
distribution is relatively insensitive to the change of yield condition, but 
the plastic zone is expected to be sensitive to the yield condition. Moreover, 
the behavior of uniformly loaded plate is more sensitive to the variation of 
yield condition than that subjected to concentrated load. 

Analysis of centrally clamped annular sandwich plates under uniform loads 
was made by French in 1964 [47], utilizing Tresca 's yield condition. The plate 
behavior has been traced from zero to the collapse load, and three phases of 
loading are considered after initial yield. Collapse pressure is given 
graphically as a function of the ratio of inner to outer radii. 

In 1964, Lackman [48] presented a method for the analysis of axi- 
symmetrically loaded circular plates. The method makes use of an analogy 
between plastic strain gradients and transverse loads. The Prandtl-Reuss 
equations together with the von Mises ' yield condition and isotropic hardening 
were employed. An example is given for a simply supported circular plate 
under uniform load. 

1.5 Experiments on Circular Plates 

Experiments other than the ones mentioned earlier to check the 
theoretical results have also been published. In 1954 Cooper and Shifrin 
[49] presented the results of nine simply supported mild steel circular plates 
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under concentric uniformly distributed load., Comparison was made with limit 
loads obtained by Hopkins and Prager [7]. Similar kinds of tests were per- 
formed by Dyrbye et al (1954) [50], Hay thornthwai te (1954) [51], and Foulkes 
et al (1955) [52] on either clamped or simply supported plates. 

Correlation of experimental evidence with collapse load using limit 
analysis turned out to be unsatisfactory. 

Hay thornthwaite and Onat 1955 [53] reported the test results of steel 
plates under reversed loading. It was observed that although the limit load 
is often of little physical significance for a monotonically increasing load, 
it becomes a measure of the minimum load carrying capacity in reversed loading. 
This question still needs further investigation. 

Lance and Onat 1962 [54] reported additional experiments on mild steel 
plates under uniform and partial load over a small circular central region. 
Comparison with collapse load in the sense of limit analysis re-affirmed the 
previous results. Etching patterns and mill-scale flaking patterns were also 
studied in this report. 
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II. ELASTIC-PLASTIC ANALYSIS OF CIRCULAR PLATES 

2 . 1 Foreword 

For the purpose of analysis, the plate is divided into a number of 
annular elements, which are further subdivided into a number of layers along 
their depths. The primary objective is to establish the stiffness matrix of 
an annular element for a typical increment of loading. The assemblage of these 
annuli can be easily accomplished using direct stiffness method of matrix 
analysis of structures. A modified analysis will also be given for a central 
circular disc element for solutions of solid plates. Loadings are applied at 
circular nodes where two consecutive elements meet. Both tributary area 
approach and consistent equivalent nodal ring load method were employed to 
convert the transverse load applied within the elements to the nodal ring 
loads. Possibility of loading as well as unloading is included in the 
analysis of each load increment. 

The following relations were established for a typical loading increment. 
Notations will be described as they first appear, and they are also collected 
at die end for reference. 

2.2 Equilibrium Equations 

By adopting the sign convention for positive quantities as shown in 
Fig. 1, the equilibrium of the increments of moments and transverse forces for 
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the case of axisymmetric loading are as follows 


d A M 
r 

dr 



- A v 


+ A Q = 0 


( 2 . 1 ) 


d A Q 1 
dr + r 


A Q = A p(r) 


( 2 . 2 ) 


Eliminating AQ 




d 2 A M 
dr 


d A M 

l 

dr 


d A M 
dr 


e 


) + A p(r) = 


0 


(2.3) 


Note that the symbol A indicates a finite increment. 

AM_^, and AM^ are increments of radial and tangential moments per unit 

length. 

AQ is increment of radial transverse shear per unit length. 


2 . 3 S train -Displacement Relations 


In accordance with small deformation theory and the assumption that 
plane section normal to reference plane before deformation remains so after, 
the strain-displacement relations in the absence of in-plane forces are 
expressed as follows 


j'A (A kfj 

V .4 - - > *3 ; e ) - - * l A 

* See for example, Timoshenko and Woinowsky-Krieger , "Theory of Plates and 
Shells,” 2nd Edition, pp. 51-53, McGraw-Hill, 1959. 
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where : 

A £ r * Ae 0 are radial and tangential strain increments, respectively. 

AK^, AKg are measures of radial and tangential deformations which are 
only functions of r. z is a coordinate distance measured positive downward 
from the reference plane. 

For axi-symmetric deformation, and neglecting shear distortion, that is 
assuming plane section normal to the reference plane before deformation remains 
normal to its deformed state, the following relations hold* 




(2.5) 


where A K , A K are now the increments of curvature of the reference surface 
r y 9 

and A w is transverse deflection increment which is taken to be positive in 
the direction of transverse load application. 

2,4 Stress-Strain Relations 

A constitutive relation for an elastic-plastic solid should include the 
state of stress and strain as well as their rates in order to account for the 
history dependence. As has been mentioned by some authors, an attempt to 
express the various properties of elastic-plastic solids by means of a single 


* Ibid 
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mathematical model can be hardly achieved. Within the framework of a 
phenomenological approach, the attempt has been made both to construct the 
general theories to be adequate for describing the phenomena occurring beyond 
the elastic limit, and to formulate it in a manner suitable for practical 
applications. Because of non-linearity and irreversibility of the deformation 
processes, even the simplest model is generally too complicated to apply. 

In order to establish constitutive relations for plastic solids three 
ingredients — yield condition, flow law, and hardening rule — must be 
defined. A number of surveys regarding the constitutive laws of plasticity 
have been published. The discussion here will be limited to those which are 
closely related to the laws chosen to demonstrate the proposed method of plate 
analysis. The general theories will not be discussed. The reader, however, 
is referred to references [55] through [62] for details of recent progress in 
this area. 

Prandtl-Reuss theory is most widely used to describe elastic-plastic 
deformation. The materials, here, are considered to be time independent and 
initially free from residual stresses. In addition, the process is considered 
to be isothermal and the deformation small. Following Prandtl-Reuss theory, 
strain tensor is assumed to consist of two components: elastic or recoverable 
and plastic or permanent. 


€ . . 

ij 




ij J — 1 > 2 ,3 


( 2 . 6 ) 


or 


d € . . = d 
iJ 


iJ 


+ d c , 


ij 


(2.7) 
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where superscripts E and P designate elastic and plastic components, respect- 
ively. Some authors, instead, have preferred to use the rate of strain 

d € . . 

^ij = — dt^ » and simila **ly rates of stress and displacement. Since, 

however, it has been already assumed that the material behavior, and con- 
sequently its constitutive relations, are time independent, the latter notation 
seems to be unnecessary. This avoids a possible confusion with viscosity, and 
its use is abandoned here. 

The elastic component of strain is related to stress through the 
generalized Hooke’s law 


d 



C ijki d T ki 


— 1 * 2,3 


which in the case of isotropic materials becomes 


( 2 . 8 ) 


d 



l+v 

E 


d t. . 
1J 


V 

E 


dT 


kk 


5. . 

i J 


(2.8b) 


Here, the repeated indices imply summation. 

In order to establish relations between the plastic component of strain, 
and the state of stress; the existence of the plastic potential, and the 
validity of normality rule at a regular point on the yield surface are 
generally assumed, that is 


e. P = d e . } 
ij ij 


df 

d T. " 

ij 


li 


= 0 


= 4<P 


( 2 . 9 ) 
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where , 

f (t € P 

v k£ 9 ki , K) = 0 

is the yield condition, and 

e - € 5., 

ij ij 3 kk ij 

is the deviatoric strain tensor. 


0 is a non-negative function which may depend on stress, stress rate, 
strain, and history of loading. 

K is a work-hardening parameter. 

The condition of loading, neutral loading, or unloading is distinguished 

df 

by whether, ^ d t . is greater, equal, or smaller than zero, respectively. 

ij 1J 

Drucker* suggested 


d<t> = G 


df.; 


d T 


ki 


ki 


( 2 . 10 ) 


where G is independent of the stress rate. 

2.4.1 Elastic-Perf ectly Plastic Solids 

Since the plastic deformation can be assumed to occur without a change 
of volume, i.e., = 0 , expression (2.7) can be re-written as follows 


d e . . = d e . . + de.. = d e. . + d<t> ^ — 

ij ij ij dt. , 

ij 


( 2 .«) 


*Drucker, D.C. , "A Definition of Stable Inelastic Material," J. Appl. 
Mech. , Vol. 26, pp. 101-161, 1959. 


Assuming the materials obey the Mises yield condition 
f = j 2 - k 2 = 0 

where 



is the second invariant of stress deviator, and 


( 2 . 12 ) 


s . . 

iJ 



kk 


5 . . 

iJ 


where k, the yield stress in simple shear is k = - — ay, where ay is yield 

•J~3 

stress in uniaxial tension. Using these definitions, 


df 

ch\ . S ij (2.13) 

iJ 

For isotropic elastic material the Hooke's law in terms of stress and 
strain deviators is expressed as 


E 1 

e. . = — s. 

iJ 2\i ij 


(2.14) 


or 



1__ 

2m* 


d s 


ij 


(2.15) 


Where, [i is the shear modulus. 

Substitutions of (2.13) and (2.15) into (2.11) yields 
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Considering that 
df = s . . d s . . = 0 


ij iJ 


multiplying (2.16) by s_, and summing, we obtain 



d0 = 


s . . d e . . 
ij iJ 


2 k 


Substituting (2.18) into (2.16) and transposing, we get 


d s. . = 2\i ( d e. . s, „ s . . de, n ) 

ij U 2 k 2 ki ij ki' 


(2.17) 


(2.18) 


(2.19) 


The hydrostatic state of stress, of the order of yield stress, does not 
produce yielding; and it is related to volumetric strain through elastic law 


T kk - <3x + ^ € kk 


( 2 . 20 ) 


or 


d T„ = (3 \ + 2|i) d €, , (2.21) 

kk kk 

where \ and \i are Lame constants. 

Multiplying (2.21) by \ 5 and adding it to (2.19), and by considering 

3 l j 

the relations between stress or strain tensor with its deviatoric components, 
we obtain 


t u - 2,ids ij 


+ \ dT 5. . 
kk ij 


% s s 
2k 2 « 13 


de 


k i 


( 2 . 22 ) 
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Note that: 


s, „ d e, „ 
ki ki 


s u WE ki - 5 ds «n V 


s de 
ki ki 


(2.23) 


Finally, 


d t. . = E. „de „ ; Idx} 

ij ijki ki .> i ^ 

9x1 


= [E] \d € ) 

9x9 9x1 


(2.24) 


where 


E iw‘* 5 ik 6 « + " 6 ij s «- ^ 


(2.25) 


\ and n can be expressed in terms of Young's modulus E and Poisson's ratio as 


. _ _ VE E 

(1+V);: (1-2V) » ** ~ 2(1+V) 


(2.26) 


Substitution of (2.26) into (2.25), gives 


E ijki l+v (8 ik 6 ji + 1-2V 6 ij 5 ki " ot 2 S ij S ki } 


(2.27) 


Generalized Plane Stress 


The general expression (2.24) can be specialized for the state of plane 
stress. For this purpose assume that t. =0, dr. =0, where i = 1,2,3 and 

lO lO 

further that de^ 3 = de 23 =0, and 


dT 33 _ E 33T& d€ r& + E 3333 d€ 33 - °* T,5 ~ 1,2 ‘ 
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Then 


-■ E 


d e 


33r& 


d€ A . 


33 E 3333 r& 


d T op = E a prs d e rs + V 33 d e 33 where a « p= 1 > 2 


Substitution of (2.28) into (2.29) gives 


d T ap E aprs d e rs 


a,P,r,s = 1,2 


and 




E aprs E 3333 E Qp33 E 33T5 
E 

3333 


For the principal directions of stress increment, where 
directions of stresses do not change during loading, the basic 
the generalized plane stress can be written down as 










d t 

11 



E mi 

E 1122 



d £ 11 

< 

d t 

22 

> 


E 2211 

E 2222 


< 

d € 22 






- 




(dir} = [E] [d e) 

2x1 2x2 2x1 


(2.28) 

(2.29) 

(2.30) 

(2.31) 

the principal 
relation for 


(2.32) 


or symbolically as 


25 


Where 


s 


= E 


22 


1111 2 2 

s, , + s„„ + 2V s s 
11 22 11 22 


E = g - E S H S 22 

1122 E 2211 E 2 2 


S 11 + S 22 + 2V S 11 S 22 


I = E 

2222 22 


'll + S 22 + 2V 'll '22 


(2.33) 


In this case, 


S 11 3 (T 11 


~ T ) 
2 ^22 J 


2 1 

S 22 " 3 ' ,T 22 " 2 T ll' 


(2.34) 


Note that the [E] matrix in (2.32) is singular. This is also true in the 
general case of (2.24). This follows from the fact that the stress increments 
d T are linearly dependent as can be seen from (2.17), since 


s . . d s. . 
ij ij 


= s . . d t. . = 0 
ij ij 


(2.35) 


Therefore, we can only specify strain increments and solve for the 
stress increments; however, the reverse is not possible. This is obvious 
for a uniaxial stress condition for a perfectly plastic solid by noting its 


stress-strain diagram. 
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The expression (2.24) is valid for loading when 


df 

3F 


ij 


d x. . = s . . d t. . = 0 
ij iJ ij 


(2.36) 


That is when the point in the stress space moves on the yield surface 
or stays on it. Stated otherwise, 


f (1) = f (2) = 0 


(2.37) 


But, if 


|£— d T.. <0 


unloading takes place, and instead of (2.24) the Hooke's law must be used, i.e., 



1 



- 

- 




d T 11 

} = 

E 

1 

V 

i 

1 d e ll 1 

d T 22 


l-v 2 

V 

1 


d e 

22 


J 






J 


(2.38) 


It is interesting to note, that for a strain increment vector, which 
is in the direction of normal to the yield surface, the stress increment 
vector vanishes. Therefore, if for a stress point located on the yield 
surface, we resolve the strain increment vector into two components — tangent 
and perpendicular to the yield surface, see Fig. 2 — it is only necessary to 
retain the tangential component in the computation of stress increment. Hence, 
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\d = [E] \de) = [E] [de T J +[E] { de N } . 

^ d t| = [E] ^de T ) 


(2.39) 


This property is preserved for the generalized plane stress case dis- 
cussed earlier. In this case both {dT^ and ^d€ } are tangent to the 

yield curve. This indicates that the transformation [E] only stretches 
T 

^d€ j with a proper scale. 

The stress-strain relations (2.32) are strictly true for infinitesimal 
increments. Here it is necessary to establish relations for finite increments. 
Thus, the infinitesimal increments are to be integrated within a small interval. 

3|C 

Prager [34] and Hodge suggested the concept of piecewise linear yield condition. 
It enables one to integrate stress-strain relations locally while preserving 
the incremental characteristics in the total. Another approach has been 
employed here. 

The equations (2.32) essentially specify an initial value problem. For 
a finite increment of strain, the stress increment should move along the yield 
curve as from A to C shown in Fig. 3. By following a tangent at A, one reaches 
an incorrect point B. However, by projecting point B back to the yield curve 
point C is located. This point C is taken as the new state of stress from 
which the next step in calculations is made. In Fig. 3 the locations of points 
B and C are greatly exaggerated and the result of example I, Fig. 14 indicates 
that for a relatively small increment of stress the distance between such 

*See footnote on pp. S. 
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points are small. Therefore the error committed in this process is not very 

significant. As the first approximation the initial state of stress in each 

increment can be used to determine the [E] matrix. Iteration procedures or 

the techniques of numerical integration such as the Runge-Kutta or Euler’s 
* 

modified method was employed to compare the results found with the first 
approxima tion . 

A difficulty arises in the transition zone, when the point initially 
located inside the yield surface reaches the surface during the increment of 
loading. In this process the associated deformation contains a purely elastic 
part. During deformation the point moves elastically within the yield surface 
until it just touches the surface. When the point reaches the yield surface 
the plastic part of the deformation occurs. If the total deformation produced 
during an increment of loading is given, equation (2.38) can be used to 
determine the corresponding stress increment. Depending on the relative dis- 
tance of the initial and final state of stress with respect to the yield 
curve, point D and E, of Fig. 3, the point closer to the yield surface, point 
E in this figure, is chosen for a radial approximation. The intersection of 
the radius vector with the yield surface, F in Fig. 3 is taken as the point 
of the initiation of yielding. Strain increment associated with DF is obtained 
using (2.38). It is subtracted from the given strain increment to find the 
strain associated with the elastic-plastic deformation. The remainder of 
the procedure follows as before. 

* See for example Collatz, L. , "The Numerical Treatment of Differential 
Equations." pp. 53-61, Springer-Verlag , 1960. 
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2.4.2 Elastic-Plastic Strain Hardening Solids 

From the observations of the behavior of hardening materials in uni- 
axial and biaxial tests, it is known that during plastic deformation the 
yield surface is continuously changing in size and shape. 

Isotropic hardening (see Fig. 4a), which at higher stresses exhibits 
uniform expansion of the initial yield surface, is the most widely used law 
to describe hardening. Thus, the yield surface for initially isotropic 

h 

materials depends on a single parameter k and may be written as 


f = g (J 2 , Jg) - * = 0 (2.40) 

where J 2> are the second and third invariants of stress deviator, and K 

is the hardening parameter which describes the strain history. Two measures 

of hardening for K are frequently used. The first approach suggests that 

the degree of hardening is a function only of the total plastic work, and is 

* 

otherwise independent of the strain path. 


K = K (W ) 
P 


w 

p 


- / 


T. . d 

ij 


ij 


(2.41) 


where is positive definite, since plastic deformation is an irreversible 
process. The second approach states that the so-called equivalent plastic 
strain increment 


d € P 



[ d € 


ij 



1/2 

] 


(2.42) 


* Hill, R. , The Mathematical Theory of Plasticity, 1 ' pp. 23-33, Oxford, 
Clarendon Press, 1950. 
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integrated over the strain path provides a measure of the plastic deformation . 
That is 


K = k (/ d e p ) . 


( 2 - 43 ) 


* 

As pointed out by Hill , the above two concepts are equivalent for the 

materials obeying von Mises yield condition Hill’s remark was later 

** 

generalized by Bland , showing that for any g to be a homogeneous polynomial 
of degree n, and have regular regimes of yield surface, the above property 
holds true, if, and only if, g is linear or quadratic function in the 
principal components of stress. Certain restrictions on the coefficients are 
imposed in the quadratic case. 

For proportional loading, if the strain ratios remain constant, (2.43) 
can be integrated 


€ P = / d i P 




1/2 



(2.44a) 


K = K (e P ) (2.44b) 

jjesjojc 

Expression (2.44b) was checked by some investigators , with the experimental 
results, where the stress ratios were not constant. In the tests on mild 
steel and annealed copper, agreements to within 5 percent are reported. 


* ibid 

**Bland, D.R. , "The Two Measures of Work-Hardening , " Proc. 9th Int. Congr. 
Appl. Mech. (Brussels, 1956), Vol. 8, pp. 45-50, 1957. 

***See Hill, R. , footnote on pp. 29. 
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The concept of isotropic hardening does not account for Bauschinger 
effect. Actually, it predicts a negative Bauschinger effect. Therefore, 
isotropic hardening would predict erroneous results in problems involving 
unloading followed by reloading along some new path. To include the 
Bauschinger effect, a hardening rule was suggested by Prager*, which assumes 
® translation of the initial yield surface (see Fig. 4b) . Prager 

employed a kinematioal model to describe this hardening rule. For this reason 
it is termed "kinematic hardening." This hardening rule can be represented 
mathematically by 


f(T ij - V = F(T ij - V 


(2.45) 


where f (t^) = 0 is the initial yield surface, and a is a tensor represent- 
ing the total translation of the center of the initial yield surface. Prager 
suggested that the yield surface be translated in the direction of the normal 
to the yield surface for any increment of strain. 


d “ij = C d € ij (2.46) 

where c is generally assumed to be constant. In such a case the process is 
called linear hardening. Shortly after Prager’s proposal it was recognized 
that the properties of preserving the shape, and of pure translation of the 


♦Prager, W. , "The Theory of Plasticity: A Survey of Recent Achievements," 
(Jame Clayton Lecture) Proc. Inst. Mech. Eng. Vol. 169, pp. 41-57, 1955 . 
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yield surface along a normal, do not generally remain invariant if the nine 
dimensional stress space is degenerated into subspaces. To resolve this incon- 
sistency Ziegler* suggested to replace (2.46) by 

da = (t - a ) dn (2.47) 

ij ij iJ 

Further suggestions such as piecewise linear yield conditions, which 
accomodate both translation and expansion of the yield surface (see Fig. 4d) , 
and the concept of yield corner stating that the yield surface changes only 
locally (see Fig. 4c), have also been advanced. Numerous tests have been 
conducted to check these theories, but no definite conclusions have been 
reached so far. 

After considering the several possibilities of a constitutive law, 
isotropic hardening was adopted in this report. The material is assumed to 
obey von Mises ’ yield condition. Then, in equation (2.40) the function g 
becomes 

g = J 2 (2.48) 

Assuming the validity of (2.44b) for non-radial loading; and sub- 
stituting (2.44b) and (2.48) into (2.40) , we obtain 

f=J -K(e P )=0 (2.49) 

Ct 

}j( 

Ziegler, H. , "A Modification of Prager's Hardening Rule," Quant. Appl. 

Math. Vol . 17, pp. 55-65, 1959. 
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The unknown function it can be determined from a simple experiment; 

namely, uniaxial tension or simple shear. The expression (2.49) has been 

- 1 *- 

sometimes called the J 2 - theory and is frequency expressed as 


o = H (e P ) 


(2.50) 


where , 


1/2 


= s , H ] ='^3jT 


IJ 1J 


(2.51) 


is the so-called effective stress. 

Now, recall (2.9) and (2.13) which gives 


d €. = d <J> s 

iJ ij 


(2.52) 


To determine d$ multiply (2.52) by - d and sum, then note (2.42) 

and (2.51) to obtain 

(d e p ) = | (d e p d e p ) = | (d<t>) 2 s. . s. . = (§ a d <P) 2 
J 3 xj ij 3 ' 


d ♦ - | 

2 - 


(2.53) 


If the curve of a - e p is given as in (2.50), then 


d e " p = 52 

H' 


(2.54) 


where H' = 


da 


d e 

into (2.53), we obtain 


— is the slope of curve. 


Whence upon substituting (2.54) 
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d <D = 3 i_ ££ 

2 H’ 5 


3 i_ 

2 H' 


c>5 


d t 




ki 


(2.55) 


Also note that 


d J„ 


ij 


37 ; 


ij 


2 - S5 

= 3 a 37“ 

ij 


(2.56) 


Substitute (2.55) and (2.56) into (2.52) to obtain 


„ _ P _ 1 
ij H’ 


Scr Scr 

37“ St 

ij 


d t 


ki 


k Z 


(2.57) 


This relation is similar to expression (2.10) proposed by Drucker. It 
can be re-cast into another form: 


d € 


3 1_ 

ij “ 2 H’ 


S ij S kl d _ 

s s kJt 

mn mn 


(2.58) 


If the data for uniaxial tension test are used to define H f , it is easily 
seen that 


1_ 

H' 



1 

E 


(2.59) 


where is the tangent modulus, and E the Young’s modulus. If the data 
from pure shear test are used, then to define H’ we have 


1_ 
H 1 


1 

3 



(2.60) 


! 
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where ^ is the tangent modulus in shear for t -T diagram, and n is the elastic 
shear modulus. 

For H ' to be invariant, the comparison of (2.59) and (2.60) leads to the 
following requirement: 


3 1 1-2V 

E t u " E (2.61) 


Expression (2.61) imposes a restriction on E t and ji t , which generally does 
not hold true for all materials. 

Now, returning to Prand tl-Reuss equation (2.7) and substituting for 

d€ ij ^ rom (2.58), and considering the uniaxial tension test as the basis, we 
obtain 


d € 


ij 


1+V 

E 


d T 


ij 


V 

E 


d t 


kk 


5. , 
ij 




S ij S ki 


S mn S mn 


d (2.62) 


or 


d € . 


ij 


S ijki d T ki 


(2.63) 


where 


S = 5 5 - X 6 6 + 3 , 1_ _ 1 , S ij S jq 

J E lk jjt E ij b ki 2 ( E E s s 

t mn mn 


i jki 


(2.64) 


Symbolically we can express (2.63) in matrix form as 


{de} = [S] (dr) 

9x1 9x9 9x1 


(2.65) 
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Specializing (2.63) for the case of plane stress in principal stresses 
and following the similar procedures for elas tic-perf ectly plastic solids, 
we obtain 


1111 S, 


Where , 

1 2 
(t - - t ) 
s = 1 / !_ Is 11 2 l 22 ; 

1111 E V E E ' 

t a * 

s -S V 1_ 1 (t 11 ~ 2 T 22 } (x 22 

1122 2211 E 1 E + E } 

t a * 

1 2 
(T - — X ) 

s = 1 . , 1_ In 22 2 ll' 

2222 EVE E ^ - 2 

t or 


1 ^ 


Inversion of (2.66) gives 


E llll E 1122 


E 2211 E 2222 
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Where 


E 


1111 


1122 


E, 


2222 



£ + 

3 


2 

S 2 

= E 

4 

2 

S 1 “ 

S 1 S 2 + S 2 2 


(1-v 2 ) 

o s 

6 + ? __ 

* + 2V 

2 

s, 

12 2 


5 4 

s 

2 

1 - S 1 

2 

2 2 

= V. 

- V 

V 

«-! 

s s 
1 2 

2 2 
s i - s i s 2 + s 2 

2211 

(l-v 2 

> 5 + ! 

S l 2 + 2V S 1 S 2 + S 2 2 



2 2 
S 1 " S 1 S 2 + S 2 



‘♦5 


2 

S 1 

= E 


2 

S 1 - 

2 

S 1 

12 2 


(1-V 2 ) 

t 3 

5 + T 

S l 2 + 2V S 1 S 2 * S 2 2 


S 1 “ S 1 S 2 + S 2 


(2.69) 


where 


£ = 


2 1 

E-E t ’ S 1 = 3 (T ll" 2 T 22^’ S 2 = 3 (T 22“ 2 T ll^ 


- 2 2 2 
CJ — T — T T +T 

11 11 22 22 


(2.70) 


It is easily seen that if E^ tends, to zero expressions (2.69) will be 
reduced to (2.33). 

Stress-strain relations (2.32) and (2.68) will be employed to describe 
the material properties of any annular layer during loading. During unloading, 
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apply relations (2.38). Since in the computation processes we are dealing 
with finite increment of loading, (2.32) and (2.68) are not strictly valid. 

As mentioned earlier, however, the errors introduced are not very significant 
for small increments. 

Although the examples which are given in this report have been based 
on (2.32) and (2.68) the developed procedures are perfectly general and any 
other constitutive relations can be employed by properly defining the [E] 
matrix. Should necessity arise, experimental data may be used directly. 


2.5 Stress Resultants 


To determine stress-resultant increments, we proceed as follows: 

Take the middle plane as the plane of reference, and divide the plate 
into a number of layers symmetrically arranged with respect to the middle 
plane (Fig. 5). Then upon assigning the material properties to each layer, 
we can formulate the expression for as 


f A f A 

|am) = 2 ) J | At) zdz = 2 ) / [E (k) ] {Ae} zdz (2.71) 


k=l k-1 


k=l k-1 


where \ Am| contains both the radial and the tangential moments, i.e., 

=, 


AM ' 

I AM 

1 

i r 


' 

AM„ 

AM 

2 

l e J 
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Upon substituting from (2.4), and integrating, we obtain 

^Am] = - [D] ^Ak] 

where, with h = o, 


[D] =| ^ [E (k) ] 


(h, 


h ll> 


k=l 


If the layers are taken of equal thickness 


(2.72) 


(2.73) 


h = k h ; h = £- 

k 1 2n 


(2.74) 


and 


n 


[D] 

2x2 


h 
12 n 


3^ [E (k) ] (3 k 2 - 3 


k+1) 


k=l 


(2.75) 


here, 


[D] = 


D D 

11 12 


D D 

21 22 


’ WhCre °12 = D 21 


(2.76) 


This equation together with (2.72) yields the following results: 


AM 

= AM = ■ 


AK - 

- D „ 

AK 

r 

1 

11 

r 

12 

e 

AM 

= = • 

- D oi 

AK - 

- D_ . 

AK 

e 

2 

21 

r 

22 

e 


(2.77) 
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2.6 Governing Differential Equation 

Substitute (2.77) into the differential equation of equilibrium (2.3), 
and define curvatures as given by (2.5). Then the governing differential 
equation becomes 


D 


11 


Aw 



III 


11 


h D 22 AW + ~3 D 22 AW ’ = A P(r) 

r r 


(2.78) 


Dividing by , and letting 


vV 2 


D 


22 

A 


IV 2 III 
Aw + — Aw 
r 


A (Aw 
r 


■=■ Aw’) 
r 


Ap(r) 

D 11 


(2.79) 


(2.80) 


where primes denote differentiation with respect to r. 

The solution of part of homogeneous (2.80) has three different forms 
depending on the value of 


(a) For _/\_ = 1 

2 2 

Aw = a, + r + a^ in r + a A r in r 
12 3 4 

(b) A= 0 

2 

Aw = a + a r + a r + a r in r 

-L ^ o 4 

(o A^i.A^o 

i+A i-A 2 

Aw = a, r + a r +a„r + a 

1 4 o 4 

Case (b) is not frequently encountered. 


(2.81) 


(2.82) 


(2.83) 
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2.7 Stiffness Matrix For an Annular Element 

In this section the relations between nodal ring forces applied at the 
edges of an annular element, and their corresponding nodal ring displacement 
are established. Here the term "nodal ring force” includes both transverse 
force and moment. Similarly, the "nodal ring displacement" implies both 
displacement and rotation. 

In order to set up the element stiffness matrix, two frequencly en- 
countered cases A= 1 andA^ 1 , A/ 0 will be considered. 

Case I A= 1 

From expression (2.81) we can obtain all nodal forces and displacements 
in terms of a. , i = 1,2, 3, 4. 

Thus, 

Increment of rotation: 

A co = A w' =2 a 2 r + a g r 1 + a 4 r (2 in r + 1) (2.84) 

Increment of radial curvature: 

A K r = A w" = 2a 2 - a 3 r 2 + a 4 (2 in r + 3) (2.85) 


Increment of tangential curvature: 


A K @ = r A w' = 2 a 2 + a 3 r" 2 


+ a 4 (2 in r + 1) 


( 2 . 86 ) 
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Substitution of (2.85) and (2.86) into (2.77) yields moments 
A M r = - 2 . 2 <D U + D 12 ) - a 3 (-D^ + D^) r' 2 

- a 4 (D n (3 + 2 in r) + D (1 + 2 in r) 

A M 8 - - 2 =2 (D n + D 12> + ° 3 ( - D u + D 12> r ' 2 

- a 4 [D X1 (1+2 in r) + (3 + 2 in r) ] 

Finally, substitution of (2.87) and (2.88) into (2.1) gives 
A Q = 4 a 4 D n r _1 


(2.87) 


( 2 . 88 ) 


(2.89) 


Next, consider the annular element shown in Fig. 6 with all quantities 
having positive sense. Nodal ring forces |A S^can be defined in terms of 
a i , i = 1,2, 3, 4, as 



[B,] \A 

4x4 4x1 


where 



A Q 


A M 


A Q 


A M 


r / 



( \ 
a. 


= 


a. 




(2.90) 


(2.91) 
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The elements of [B ] are given in Appendix A. Note that because of 
the change of sign convention from that of strength of material shown in Fig. 
1, to the one shown in Fig. 6, the following modification should be made: 

AM (r .) = - A M J 
r J r 

(2.92) 

A Q (rj = - A Q J 


Similarly, nodal ring displacements can be expressed in terms of 
a i , i = 1,2, 3,4, as 


{Av} = [BJ. M 

4x1 4x4 4x1 


where, ^aj 



is defined as before in (2,91), and 

i \ 


A w 


A CD 


A w" 


\ , A cd 1 ’ J = (A w ' ) 


1 > J 


A CD J 


(2.93) 


(2.94) 


The matrix [B^] is given in Appendix A. To evaluate the stiffness matrix, 
solve (2.93) for ^ aj , and substitute in (2.90). Thus 


M = t ^]" 1 (A vj 

4x1 4x4 4x1 


(2.95) 
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, i " 1 

} A s\ = [B g ] [B y ] \A v\ = [k] $AvJ 
4x1 4x4 4x4 4x1 


where [k] is the element stiffness matrix, and 

-1 

[k] = [B ] [B ] 

s V 

4x4 4x4 4x4 


(2.96) 


(2.97) 


Case II , A^o. i 

In this case nodal forces and displacements can be derived from 
expression (2.83), similar to that of Case I: 

Increment of rotation: 


A oo = A w ' 


B 1 (1 +A) 



(i -A) 


-A 

r +2 a g r 


(2.98) 


Increments of curvatures 

A K r = A w" = a x (1 + A> - a 2 A(l"A) r"^'" 1 + 2a 3 

A K 'Q = r A w ’ = a 1 O +A) r^" 1 + a 2 (1-A) r" - ^" 1 + 2a g 

(2.99) 


Moments -- Substitute (2.99) into (2.77) to get 

AM- = -a, (1+A) (AD n ,+D n0 ) 7 /^ -a 0 (1-A) (-Al> .,+D ) r 


11 12 ' 


11 12 


- 2a 3 (D 11 +D 12 ) 


( 2 . 100 ) 
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and 


A M = -a., (1+A) 


<d 22 *AV 


A-i 


a 2 (1-A) (-D 22 -«Ad i2 ) 


-A-i 


-2 a 


3 


(D„ „+D ) 

12 22 


( 2 . 101 ) 


Utilizing the equilibrium equation (2.1), upon substitution into it of 
(2.100) and (2.101), we obtain 


A Q = 2 a 3 (D n -D 22 ) r 1 (2.102) 

Following the same steps as shown in (2.90) through (2.96), we obtain 
the stiffness matrix for Case II. Formally, the equations are the same as 
before, 

|A s] = [k] ^ A v] 

4x1 4x4 4x1 

and 

[k] = [B s ] [B y ] 

4x4 4x4 4x4 


(2.96) 


(2.97) 


However, the matrices [B ] and [B ] for this case are different as can be 

S V 

seen in Appendix B. 

2 . 8 Stiffness Matrix for a Disc Element 

The condition of axial symmetry, and the requirement that the load is 
applied only at nodal rings require that A Q for a disc element should vanish 



throughout (see Fig. 7). For such an element two cases of different -A-'s 
must again be recognized. 


Case I W= 1 


(1) Since A w must be finite, a 3 = 0 


(2) Since AQ= 0, a„ = 0 

4 


Therefore, as before, 


. 2 
A w = a + a r 
1 ^ 


A a) = A w 1 = 2 r 
2 


2 

r . 

J 


2r . 
J 


2 

r 

J 

2r . 

J 


-1 

[B ] 

X> 


Also 





(2.103) 

(2.104) 


(2.105) 


(2.106) 


(2.107) 


A M - - 2 a. (D 11+ D ) = A M 


(2.108) 
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A Q = 0 


Hence , 


( 2 . 109 ) 




0 




2 < 1) 11 +D 12> 



( 2 . 110 ) 


That is 


fV 

2x2 



0 



2<D 11* D 12 ) 


(2.111)o- 


Theref ore 


[k] = 

= [ v 

[B ] 
L 0 J 

2x2 

2x2 

2x2 


(D +D ) r 
11 12' j 




(2.112K\. 


The Case II, 1,A ^ 0 is not encountered in this problem since 

the ratio of radial to tangential moment increment remains constant (i.e. , 
unity). Therefore, this case will not be discussed here. 

Note that stiffness matrices developed for annular and disc elements 
are singular. This condition is due to the fact that the elements of vector 
must satisfy the condition of equilibrium, and are therefore not 
linearly independent. 


\ 
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2.9 Plate Stiffness Matrix 

The element stiffness matrices developed in Articles 2.7 and 2.8 can be 
easily assembled to obtain the plate stiffness matrix, by considering the 
conditions of equilibrium, and compatibility at nodal rings. 

Consider the n-th element, and partition its stiffness matrix as shown 


'as 1 ' 

n 


J 

k. n 

11 

1 

1 

1 

1 

1 

k. n 

iy 
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y 

A S'-* 
n 

A -tr 1 


k 11 
. J1 

1 

1 

1 

t 

A A 

k n 

jj 



4x1 4x4 


A v 


A v J 
n 

4x1 


( 2 . 111 ) 


At a nodal ring, where the n-th and the (n+l)-th elements meet, the 
equilibrium equations are 


(ah ] 

, = <Ja S J 

> + < 

A S 1 1 

"J 

1 n 1 


n+1 

2x1 

1 *■ J 

2x1 


L J 

2x1 


( 2 . 112 ) 


Compatibility of deflection and slope at the same nodal ring require 


that 


11 

u 

<1 

A v j i 



Av j 

n J 1 

n 


l H 

2x1 

^ > 
2x1 


2x1 * 


(2.113) 
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where 


( a M 

H) 


is plate (external) nodal ring force. 


is plate nodal ring displacement 
From (2.111) it is seen that 


l AS n J } = [k ji ] l A \} + [k jJ ] { Av n} 

Similarly, for (n+l)-th element 

1 AS Li} = [k "t^ { A Vl} + t k ”j 1 ] {AvU (2.115) 

Substitution of (2.114) and (2.115) into (2.112), and consideration of 
(2.113), yields 


i AR n} "t k jJ { A r n-l} * '■ k jj + tO l A V(} 

2x1 2x2 2x1 2x2 2x1 2x2 2x1 

(2.116) 


repeating the above procedure for all the elements, we obtain stiffness matrix 
for the entire plate, which is a tri-diagonalized matrix of the assemblage. 
This matrix can be schematically represented as follows: 


(a) [K] 

2Nx2N 



(b) |a r} = [K] {a r} 
2Nxl 2Nx2N 2Nxl 
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where N is the number of elements. Note that the stiffness matrix [K] 
established here is singular. It is necessary, at this point, to impose the 
boundary conditions, and to make the proper deletions in [K] in order to make 
this matrix non-singular. Then, the remaining set of equations (2.117 b) can 


be solved for 


A r 


H 


, moments, curvatures, 


Having obtained 

and stresses can be found using appropriate expressions established before. 

Since, in this report, nodal ring forces per unit length were taken, 
the stiffness matrices developed here are not symmetric. It is easy to show, 
however, that they satisfy Betti’s law. 


2.10 Transformation of Distributed Transverse Loads on the Element to Nodal 

Ring Forces 

Both tributary area and consistent equivalent nodal ring forces were 
used in the solution of problems. 

The tributary area approach simply concentrates the distributed loads 
at nodal ring junctures. The distributed loads extending half-way between the 
neighboring elements are included in the nodal load. 

Consistent equivalent nodal ring forces are derived by equating the 
virtual work of transverse loads on the element through the given displacement 
pattern to virtual work of their corresponding nodal ring forces. This is 
briefly explained below. 

Expressions (2.81) or (2.83) can be written as follows: 

T 

A w = | a j { <t, m ( r )j .01=1,2,3,4 (2.118) 

1x4 4x1 

♦Archer, J.S., "Consistent Matrix Formulations for Structural Analysis Using 
Finite-Element Techniques." AIAA Jour. Vol. 3, No. 10, pp. 1910-1918, Oct., 1965. 
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Recalling (2.93), we also have 
l A ’} ' tB v> I ■'} 


(2.93) 


Note that a bar over the symbols in the above equations is used to 
designate the quantities associated with the virtual displacements. 

The work of a distributed transverse load A p(r) on the element through 
the virtual displacement Aw is 


«, - /: 


A p(r) A w (r) 2itrdr = 2% 


t . 
1 


W’{ 


A p* 


m 


(2.119) 


where 


{ Ap *m} - J r> A V (r) 


m = 1 ,2 ,3 ,4 


( 2 . 120 ) 


On the other hand, the work of equivalent (consistent) nodal ring forces 
{APj through the virtual displacement |a vj is given 


as 


) T J 

[ 2 “ r i Ap il 

= 2it | 

.1 T T , 

" I [B V ] 1 

[ r i AP i 1 

(2.121) 

J 1 

L 2n I 

LT 

A P f 




J 

j 1 



Equating (2.119) and (2.121), and considering that the generalized 
coordinates |a|. are linearly independent, we get 


r j Ap j 


t B v T i 


-1 


{ Ap :i 


( 2 . 122 ) 


4x1 


4x4 


4x1 
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where Pj are the forces per unit length. 

Having established (2,122), the nodal ring forces are obtained from the 
following relation 


i A r „} = 

2x1 



i Ap ”) 

2x1 


+ {a p- 1 } 

2x1 


(2.123) 


where , 


r 

n 


n n+1 

r . = r . 

J i 


In (2.120) any distribution of A p(r) over the element can be specified. 
It is also possible to develop expressions for linear or parabolic distribution 
of transverse load acting on the element, and any distribution can then be 
approximated by a linear or a parabolic curve. The linear approximation was 
used in this investigation, and the appropriate expressions for j^A p | are 
given in Appendix C. 

The above procedure can be used to transform transverse ring load on an 
element to nodal ring force as well. Perhaps, however, it is preferrable in 
such cases to take the applied ring load as the nodal ring load. 

A comparison of the results using the two approaches discussed above 
of transforming the transverse distributed forces to nodal ring forces shows 
that the difference is not very significant for the small size element. 
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2.11 General Remarks 

The procedure employed in the analysis of plates described in this 
report can be outlined as follows: The plate is assumed to be initially free 

from residual stresses. The first increment of loading is applied and the 
magnitude of the load increased so that yielding just begins at some points 
within the plate. The loadings is then continued in small load increments. 

For each increment, after the Eqns. (2.117b) have been solved for | A rj , 
curvatures, strain increments, and stress increments are determined successively. 
The new state of stress is found, and [E]° matrices are computed. For the 
first approximation [E]° are taken to be the material properties for the next 
increment. In a more refined procedure, Euler's modified method is utilized. 
That is, the [E] is utilized to calculate r| for the same increment. 

The new [E] matrices are established, and the average of these two is taken 
to be representative of the material properties for the increment in questions. 

Depending on the magnitude of the load increment, this modification was 
relatively of more significance for perfectly plastic materials than those 
exhibiting hardening. It was observed that the first approximation was satis- 
factory for a relatively small loading increment and hardening materials. 

The influence of the magnitude of the loading increments was studied, and some 
results are shown in the next chapter. 

A computational difficulty arises at the instant of initiation of the 
plastif ication of an element. If the radial and tangential stress increments 
differ slightly, the material exhibits a slight anisotropy after exceeding the 
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elastic limit. The corresponding changes from unity to a value very 

close to unity. Therefore, the matrices [B^] 1 and [k] become ill-conditioned. 
To overcome this difficulty two approaches have been used. The first makes 
use of double precision computation procedure , which retains 16 significant 
digits in the analysis. The second consisted of expanding the elements of 
[B ] 1 and [k] matrices in terms of (1-/V), clearing fractions, and treating 
the ill-conditioning factors separately. The expanded forms of elements of 
these matrices are given in Appendix D. 
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III. APPLICATIONS 

This chapter consists of two parts. In the first part several examples 
are presented; in the second, a description of computer programs is given. 

3 . 1 Examples 

Eleven examples have been prepared to illustrate the main features, and 
to test the accuracy of the method of solution developed for circular plates 
in Chapter II. 

Except for Example 11, the same size plate is used throughout to 
facilitate comparisons for different boundary conditions and also for different 
material properties. The total magnitude of load applied to the plate was 
controlled by the deflection of the plate. It is believed that in using the 
small deformation theory of plates a total central deflection of up to 0.3 
the plate thickness does not induce appreciable membrane action. The developed 
solution is applicable to relatively thick circular plates. 
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Example 1 

The behavior of the simply supported plate shown in Fig. 8 under uniformly 
distributed load is analyzed. 

I d = 20. 0 in. I 



7 

The material properties are: E = 10 psi, V = 0.24. The yield stress 

in uniaxial tension a ' = 16 ksi. The material is elastic perfectly plastic. 

The number of elements and layers is 20 and 40, respectively. Load 
increments of 5 and 4 psi are used. 

The results are plotted in Fig. 9 to 14 . Comparison with a solution due 
to V.V. Sokolovskii who used the deformation law of plasticity is given in Fig. 
15. 

Example 2 

The same plate as in example 1 but with clamped outer support is 
studied . 

The numbers of elements and layers are 20 and 40, respectively. Load 
increments of 10 psi are used. 


Results are plotted in Fig. 16 to 19. 


RESIDUAL DEFLECTION AFT 



DISTRIBUTIONS OF DEFLECTION 
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FIG. 14 STRESS PATHS 
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FIG. 15 COMPARISON WITH SOKOLOVSKII'S DEFORMATION SOLUTION 
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FIG 17 DISTRIBUTION OF RADIAL BENDING MOMENT 
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FIG 19 ELASTIC- PLASTIC BOUNDARIES 
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Example 3 

The purpose of this example is to study the effect of history of external 
loading on the behavior of the simply supported elas tic-perf ectly plastic plate. 

The plate is exactly the same as in Example 1 except that V = 0.3. 

First the plate is loaded as shown in Fig. 20-a to obtain the maximum 
possible elastic deformation. Thereafter, two different schemes of loading are 
used: In the loading sequence 1, triangular increments of load are added to 

reach the final load of Fig. 20-d. In the loading sequence 2, load increments 
are added from r = 3.75 in. to the outer edge of the plate until the final load 
of Fig. 20-d is achieved. 


Results are plotted in Figs. 21 and 22. 







0 0.2 0.4 r 0.6 0.8 1.0 

" cf 

FIG 22 FINAL DISTRIBUTIONS OF RADIAL MOMENT DUE 
TO DIFFERENT LOADING SEQUENCES 
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Example 4 

The purpose of this example is to study the effect of load history on 
the behavior of clamped elastic-perfectly plastic plate. 

The plate is exactly the same as in Example 1 except that V = 0.3 

After the first step of loading (Fig. 23-a) two different loading 
schemes are used to reach finally the load in Fig. 23-d. In loading sequence 
1, typical load increments are as in Fig. 23-b. In the loading sequence 2 
load increments are added from r = 4.75 in. to the boundary of the plate 
(Fig. 23-c) . 

The results are plotted in Fig. 24 to 26. 
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FIG. 25 DISTRIBUTION OF MOMENTS DUE TO 
DIFFERENT LOADING SEQUENCES 
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Example 5 

The purpose of this example is to show the effect of the magnitude of 
load increments on the moments and deflections. The numerical experiments are 
on 0.75 x 16 in. Simply supported elastic-perfectly plastic plates with E = 
10.6 psi and or^ = 16,000 psi. The numbers of elements and layers are 20 and 
40, respectively. It should be noted that these results were obtained before 
including the Euler modification to take into account the variation of material 
properties within a loading step. 
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FIG. 27 VARIATION OF DEFLECTION AT r = 0 USING 

DIFFERENT LOAD INCREMENTS 
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Example 6 

A 1x20 in. simply supported plate is subjected to uniformly distributed 
load. The material property of the plate is shown in Fig. 30. V = 0.33. 
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Fig. 30 

Uniaxial Stress-Strain 
diagram 

18 elements are used as indicated in Fig. 29. The thickness is 
divided into 40 layers. Load increments of 15 and 10 psi are used. 


The results are plotted in Figs. 31 to 35. 















88 


Example 7 

A clamped plate of the same size and material property as in Example 
6 is subjected to uniformly distributed load. 
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Fig. 36 


There are 40 layers; and 18 elements which are distributed as shown 
in Fig. 36. Load increments of 25 and 20 psi are used. 

Results are plotted in Figs. 37 to 40. 
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FIG. 39 DISTRIBUTIONS 0 
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Example 8 

A simply supported plate exactly the same as that in Example 6 except 
with the uniaxial stress-strain diagram shown in Fig. 46 is subjected to 
uniformly distributed load. 

The E^-a diagram of the stress-strain curve is shown in Fig. 41. The 
reason for using such a stress-strain curve is that it has an initial linear 
elastic part and also it is quite easy to modify the two other linear parts 
of the E^_-a diagram to obtain different shapes of the stress strain diagram. 

The results are plotted in Figs. 42 to 46. 



Fig. 41 

Example 9 

A clamped plate of same geometrical and material properties as that 
in Example 8 is subjected to uniformly distributed load. 


Results are plotted in Fig. 47 to 50. 
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FIG 49 DISTRIBUTIONS OF TANGENTIAL MOMENT 
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Example 10 

The purpose of this example is to compare the results of the elastic- 
perfectly plastic plate solution of Example 1 with a solution for a plate 
with hardening material having uniaxial stress-strain diagram of Fig. 51. 

Such a diagram closely approximates that of elastic-perf ectly plastic material. 

The number of elements, layers and magnitude of load increments are 
the same as in Example 1. 

Results are plotted in Fig. 52 to 55. 


<7(Ksi) 



Fig. 51 






= 4.75 
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Example 11 

The simply supported plate shown in Fig. 56 is subjected to uniformly 
distributed load. 



Fig. 56 


The material properties are tabulated in Table 1 (from Ref. 48 ) 

The plate is divided into 16 elements and 40 layers. Load increments 
of 40, 20, and 10 psi are used. 

- The results are compared with V. Lackman's solution in Fig. 57 and 58 

Table 1 



LACKMANS ELASTIC PLASTIC SOLUTION 



FIG 57 DEFLECTION OF CIRCULAR 
PLATE AT p = 300 psi 
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3-2. Description of Programs 

Two computer programs have been prepared, in Fortran IV language, for 
the solution of clamped and simply-supported circular plates subject to axi- 
symmetric loading. One is for elastic-perfectly plastic material and the 
other for hardening material. The general procedure of solution is basically 
the same in both cases, although there are differences in details due to 
limitations in numerical accuracy of solution and inherent difference in 
material properties. A discussion on the computation problems and consideration 
of the variation of material property within each step of loading is given in 
the appendices E and F. The concise flow diagrams of the main program routines 
are presented in Appendix G. A brief description of the program follows. 

(1-a) Elastic Perfectly Plastic Analysis of Circular Plates (EPPACP) 

This routine is for matrix solution of elastic perfectly plastic 
circular plate. Presently a 30 element 60 layer plate can be handled but by 
decreasing the number of layers up to 35 elements or more can be used. Nodal 
ring loads in this routine are tributary. 

(1-b) Formation of Average Material Property Functions (AMA FUN) 

The functions of this subroutine is to form the average material property 
function within each step of loading. 

( 1— c) Formation of Initial Material Property Functions (IMAFUN ) 

In this subroutine the initial material properties for the next step 
of external load are calculated. 
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(2-a) Elastic Plastic Analysis of Circular Plates (EPACP) 

This routine is for matrix solution of circular plates of hardening 
material. It is a double precision program. The number of elements and 
layers which is presently handled by this routine is 20 and 40, respectively. 
Other combinations of elements and layers can be used. 

(2-b) Equivalent Nodal Ring Load (ENRIL) 

This double precision subroutine is used to find the consistent nodal 
ring loads assuming linear variation of loads between the nodal rings. If 
it is desired to use the program for isolated concentrated ring loads, this 
subroutine should, be deleted from (EPACP). 

(2-c) Formation of Material Property Functions (MATFUN) 

In this subroutine both the average material property functions in 
each step of load and the initial functions for the next increment of load 


are computed. This subroutine is in double precision. 
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IV. CONCLUSIONS 


From the study of the examples in Chapter III as well as other observations, 
the following conclusions were reached on the calculated behavior of elastic- 
plastic circular plates. 

The comparison of elastic and inelastic moments in Example 1 indicates 
some re-distribution of moments in the inelastic state. This can also be seen 
from the distribution of residual moments in all of the examples. The re- 
distribution patterns are in general similar for elas tic-perf ectly plastic 
and hardening material. 

The stress paths in Fig. 14 indicate that the internal loading does not 
progress proportionally once the elastic limit is exceeded. Therefore, the 
deformation theory is not strictly applicable. However, results of Examples 
3 and 4, which show the influence of different loading histories, indicate 
that this effect does not contribute significantly to the general behavior 
of the plate. 

In Fig. 14 the small variation of the stress paths inside the yield 
surface justifies the application of a correction based on a radial approxi- 
mation explained in Art. 2.4.1. It can also be observed that when the increments 
of external load are small, the distance between consecutive points of stress 
paths on the yield surface is small enough to justify the tangential approxi- 
mation for neutral loading explained in Art. 2.4.1. 

The comparison between the central deflections of the simply supported 
elastic-perf ectly plastic plate as calculated in the present incremental 
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appraoch and the deformation solution due to V.V. Sokolovsky [2], Fig. 15, 
indicates the close agreement of the two solutions for loads not far above 
elastic limit. However, the incremental solution predicts greater deflections 
near collapse load. 

The results in Examples 3 and 4 indicate that moments and deflections 
in the circular plate are affected very little by different loading histories 
for the two different sequences of loading. However, the inelastic zones are 
appreciably different near the boundary of the clamped plate. 

From the results of Example 5 we conclude that as the magnitude of the 
load increments in calculations decreases, the deflections of the elastic- 
perfectly plastic plate increases. Figure 28 shows two groups of curves for 
the convergence of central moments. The full line was obtained by integrating 
the stress distribution across the depth of the plate. These moments, in 
general, satisfy the yield condition. The dashed line was determined by 
utilizing the stiffness matrix. The moments obtained in this manner satisfy 
equilibrium, but may not satisfy the yield conditions. However, as the 
magnitude of the load increments decreases they approach each other. These 
results establish bounds on the moments. It appears that the average of the 
two solutions is close to the correct one and can be taken as such. The fact 
that the variations in moment and deflection in Fig. 27 and 28 are almost 
linear can be used to advantage to predict the results when the size of load 
increments approaches zero. The dependence of the moments and deflections on 
the magnitude of load increments in hardening material varies with the hardening 
rate and is relatively less significant. 
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A comparison of the input data stress-strain diagrams, a , e, for the 
hardening materials with the one computed from the calculated state of strain 
in the plate is presented in Figs. 35, 46, 55 for Examples 6, 8, 10. A good 
agreement between the two curves is evident which indicates good accuracy and 
convergence of the developed solution for the hardening material. 

The results in Example 10 show that the solution for a hardening material 
with very low hardening rate (Fig. 55) agrees well with that for a corresponding 
elastic— perfectly plastic material. This is as could be expected intuitively. 
Part of the small differences is due to the different level of accuracies in 
the computer programs used for the two materials — double precision for 
hardening material and single precision for elastic-perf ectly plastic material. 

The comparison of the solution for hardening material presented in 
Example 11 with the one published by V. Lackman [48] indicates that for the 
uniformly distributed load of 300 psi the deflections and stress distributions 
along the depth of the plate are practically the same. 

The variation of stresses along the depth of plates indicates that the 
behavior of plates is not highly sensitive to the number of layers used in 
calculations within a reasonable range. For example, for hardening material 
of the type used in Example 11 results using 40 layers agreed well with those 


of 20. 
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Appendix A 
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Elements of [B 1 Matrix 
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Appendix B 


Elements of [B ] Matrix 
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Elements of [ b^J M a trix 
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Appendix C 

Transformation of distributed load to normal ring forces. 
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Expanded form of [B] = [By] 1 matrix, case A* i.o 
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Appendix E 

Consideration of Numerical Errors 

The numerical inaccuracies are due to two main causes both of which 

exist in matrix [B^], see Chapter II. It can be observed that the 1st and 3rd; 

and the 2nd and 4th rows of this matrix respectively are of the same form. 

For a ringe element whose inside and outside radii are of nearly the same 

magnitude, the matrix becomes ill conditioned for inversion. The other and 

more important cause of inaccuracy is due to parameter TV as it approaches 

unity A^ 1. Then some of the terms in the inversion of the [B 1 matrix take 

v 

the form of (-) which lead to computational errors. 

In general, these two causes of numerical errors do not occur in the 
same element. A graph showing the plot of j 1_ A | vs. T ± /r j for the elements, 
and the probable zone of computational error is shown in Fig. 60. This is a 
result of many tests on simply supported and clamped plates using the single 
precision computer program. 

Two measures have been taken to prevent these errors: In the program 

for the elastic-perfectly plastic material, the elements of [B y ] -1 and [k] 
have been expanded in series form as presented in Appendix D . in 
calculations enough terms of this expansion have been taken to assure 
accuracy of results. This eliminates the inaccuracy due to _/\_ but not due 
to r Vr ^ as it approaches unity. The method used to minimize computational 
errors for hardening material consisted of using double precision in the 
computer program. The results have been checked by testing whether [B'^j [b j _1 
= [''I'-]- This procedure yielded satisfactory results. 
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Appendix F 

Consideration of Variation of Material Properties Function Within A Step of 
External Load ' ~~ 

Considering the plate load and the state of stress as the independent 
and dependent variables respectively, we can write 

d t = F (t) d p (1) 

where F(t) represents a function that transforms the external loads into 
internal stresses and is expressed as a function of the state of stress. 
Equation (1) is solved numerically by replacing d p and dx by finite 
increments A p and At 

At = F(t) Ap (2) 

To solve Equation (2) with reasonable accuracy, Euler’s modified 

* 3 

method is used where the order of error is O(Ap) . 

Figure (59) is the flow diagram of the substeps taken in the computer 

program to complete the (k+l)-th step of calculation represented by Equation 

(2) for the hardening material. The procedure is basically the same for the 

elastic perfectly plastic material. The notations are the same as in 

Chapter II. 

Whenever unloading from a plastic state takes place in a layer the new 

elastic properties are utilized and this modification is not used. 

*Levy , H, ; and Baggott, E.A., " Numerical Solutions of Differential Equations , 
pp. 92-96, Dover 1950. " 


( Initial properties, known) [Ej 


TEMPORARY 
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LOADING 


ACTUAL 
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LOADING 



FIG. 59 
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CONCISE FLOW CHARTS OF COMPUTER PROGRAMS 

ELASTIC- PLASTIC ANALYSIS OF A CIRCULAR PLATE WITH AX IS -SYMMETRIC LOADING 


(EPACP J 


I* 



* THE NUMBERS REFER TO FORTRAN NUMBERS 
STATEMENTS IN THE PROGRAM . 

** IN EPPACP, THE LOAD DISTRIBUTION ON NODAL 
RINGS IS TRIBUTARY HENCE ENR1L IS 
NOT USED 

*** IN EPPACAP THERE ARE TWO SUBROUTINES 
FOR CONSTRUCTION OF MATERIAL PROPERTY 
FUNCTIONS, AMAFUN AND IMAFUN 
WHENEVER THE INCREMENT OF EXTERNAL 
LOADING IS TEMPORARY AMAFUN IS USED 
AND IF NOT IMAFUN. 


END 
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SUBROUTINE FOR THE FORMATION OF MATERIAL PROPERTY FUNCTIONS 
(FOR HARDENING MATERIALS ) 


SUBROUTINE 
MATFUN , 


CALCULATE THE DISTRIBUTION OF CHANGE 
OF CURVATURE OF THE PLATE 


If TEMPORARY INCREMENT OF EXTERNAL LOAD 


EVALUATE THE TANGENTIAL MOMENT 
AT EACH NODAL RING 


CALCULATE THE STRESS INCREMENT 
OF THE LAYERS IN EACH ELEMENT 


If TEMPORARY INCREMENTS 
OF EXTERNAL LOAD 


X- .. . 
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05 

CALCULATE THE INITIAL MATERIAL 

CALCULATE THE AVERAGE VALUE 

PROPERTIES FOR THE NEXT STEP 

•OF MATERIAL PROPERTIES IN 

OF EXTERNAL LOAD 

THIS STEP OF EXTERNAL LOAD 



1. TANGENTIAL MOMENTS 

2. CURVATURES 

3. RADIAL AND TANGENTIALS 

STRESSES OF EACH LAYER 

4. EFFECTIVE STRESS IN EACH LAYER 

5. EFFECTIVE STRAIN 


RETURN 
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SUBROUTINE FOR THE CALCULATION OF AVERAGE MATERIAL PROPERTY FUNCTIONS 
{FOR ELASTIC- PERFECTLY PLASTIC MATERIAL) 
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SUBROUTINE FOR THE CALCULATION OF INITIAL MATERIAL PROPERTY FUNCTIONS 
OF THE NEXT INCREMENT OF LOAD t FOR ELASTIC PERFECTLY PLASTIC MATERIAL) 




>0 


C 


TEMPORARY YIELD FUNCTION 
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> 


<0 


START OF PLASTIC FLOW; 
OR NEUTRAL LOADING 


ACTUAL TOTAL STRESS TEMPORARY 
TOTAL STRESS; FIELD FUNCTION - 0 


CALCULATE PLASTIC MATERIAL PROPERTY 
FUNCTIONS FOR NEXT INCREMENT OF 
EXTERNAL LOAD 
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PERFORM TANGENTIAL APPROXIMATION 
TO LOCATE THE FINAL POSITION ON 
THE YIELD SURFACE 


CALCULATE INCREMENT OF STRESS , 
TOTAL STRESS AND PLASTIC MATERIAL 
PROPERTY FUNCTIONS. 

SET YIELD FUNCTION = 0 


CALCULATE 


TANGENTIAL 

AND RADIAL 

MOMENTS 


PRINT 


1. RADIAL AND TANGENTIAL MOMENTS 

2. CURVATURES 

3. RADIAL AND TANGENTIAL STRESSES, 


AND YIELD FUNCTION FOR 
EACH LAYER 
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ELASTIC LOADING; OR UNLOADING 
FROM PLASTIC STATE 


ACTUAL TOTAL STRESS * 
TEMPORARY STRESS; 

YIELD FUNCTION ^TEMPORARY 
YIELD FUNCTION 


CALCULATE ELASTIC MATERIAL 
PROPERTY CONSTANTS FOR 
NEXT INCREMENT OF 
EXTERNAL LOAD. 
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E 

ijki 

f 


g 

G 

h 

h ] 

H 


J 2’ J 3 


M , M 
r e 

N 

p(r) 

Q 

r,0,z 


Nomenclature 

constants of integration, see (2.81) to (2.83) 
as defined in (2.46) 

elastic moduli 
deviatoric strain tensor 
Young’s modulus 

tangent modulus in uniaxial test 
elastic-plastic modulus 
yield function or loading function 
a function as defined in (2.40) 
non-negative function defined in (2.10) 
plate thickness 

the distance of k-th layer from the reference plane 

function of equivalent plastic strain defined in (2.50) 

second and third invariants of deviatoric stress tensor 

yield stress in simple shear 

radial and tangential curvatures 

radial and tangential moments per unit length 

number of elements 

transverse load per unit area 

radial transverse shear per unit length 

coordinate axes of plate, z is measured positive downward 
from the reference plane 
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s . . 
ij 

S ijki 

w 


W 

P 

a . 

ij 


P 

5 

A 



ij 

n 

K 

\ 

A= & 
V D 11 


deviatoric stress tensor 
elastic-plastic compliance 

vertical displacement of the middle plane, measured positive 
downward 

total plastic work per unit volume 

a tensor representing the total translation of the center of 
initial yield surface, see (2.45) 

defined in (2.69) 

defined in Appendix (D) 

as a prefix desiginates finite increment 
Kronecker delta 

equivalent plastic strain increment defined in (2.42) 
radial and tangential strains, respectively 
strain tensor 

parameter defined in (2.47) 
work hardening parameter 
Lame constant 

square root of the ratio of tangential to radial flexural 
rigidity of plate 




V 

6 

o 

y 

a 


shear modulus 

tangent modulus in simple shear 

Poisson’s ratio 

defined in (2.69) 

yield stress in uniaxial tension 

effective stress defined in (2.51) 
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T 


ij 


4 > 


03 


u 

t ] 

w 

4x1 

£ B s ] 

4x4 

[B ] 
v 

4x4 

[D] 

2x2 

[E] 

[E (k) ] 

2x2 

[k] 

4x4 

w 

2Nxl 

{»} 

2Nxl 

{ S 1 

4x1 


stress tensor 

non-negative function defined in (2.9) 
slope 


Vectors and Matrices 

column vector 
square matrix 

arbitrary constants, see (2.81) to (2.83) 
as defined in (2.90) 

displacements transformation matrix, defined in (2.93) 

flexured rigidity of plate, defined in (2.72) 
elastic-plastic moduli 

elastic-plastic moduli associated with k-th layer 

element stiffness matrix 

generalized nodal ring load, see (2.120) 
displacements and rotations at nodal rings 

nodal ring forces per unit length of plate 
nodal ring forces per unit length for an element 
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[S] 

9x9 

(v) 

4x1 

w 

H 

Subscripts 
a,P,r,5 
i- > J > ^ > & 


elastic-plastic compliance 

nodal ring displacement for an element see (2.94) 

strain increment vector (2.24) 
stress increment vector (2.24) 

range 1 , 2 
range 1, 2, 3 


